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Abstract 

66 '. 

, The periodic standing wave method studies circular orbits of compact objects coupled to helically 

. symmetric standing wave gravitational fields. From this solution an approximation is extracted for 

' the strong field, slowly inspiralling motion of binary black holes and binary neutron stars. Previous 

work on this project has developed a method using a few multipoles of specially adapted coordi- 
nates well suited both to the radiation and the source regions. This method had previously been 
applied to linear and nonlinear scalar field models, to linearized gravity, and to a post-Minkowski 
approximation. Here we present the culmination of this approach: the application of the method in 
^■f^ , full general relativity. The fundamental equations had previously been developed and the challenge 

' presented by this step is primarily a computational one which was approached with an innovative 

technique. The numerical results of these computations are compared with the corresponding results 
I I , from linearized and post-Minkowksi computations. 
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' I. BACKGROUND AND INTRODUCTION 

^> ' Due to successes in numerical relativity, there is now a good understanding of many features of the gravi- 
, tational waves from the inspiral and merger of comparable mass black holes [if [2, H, 0, 0, 0, H, H, [13, [Hlj [l2j 
^ ■ m, [11, [H, [H, m [3 There are, however, areas that are still best investigated with an approximation method 
for the slow inspiral epoch. One area is the evolution and gravitational waveform during the slow inspiral 
(but strong field) epoch, since there are too many orbits for numerical relativity to be feasible. Another area 
is the case of mass ratios too large for numerical relativity, but not large enough for particle perturbation 
methods. A slow inspiral approximation method is also useful as a starting point for a better understanding 
of the radiation reaction fields acting on an inspiralling hole. 

The Periodic Standing Wave (PSW) project is meant to provide such an approximation. This method 
seeks a numerical solution for a pair of sources (black holes, neutron stars) in nondecaying circular orbits 
with gravitational fields that are rigidly rotating, that is, fields that are helically symmetric. Because the 
; ^ ■ universality of gravitation will not permit outgoing waves and nondecaying orbits, the solution to be computed 
I is that for standing waves. An approximation for slowly decaying orbits with outgoing radiation is then 
extracted from that numerical solution. 

This work has progressed through several stages. In the first stage[T9l [20l [2ll [2^ . [2^, a nonlinear scalar 
fields model was investigated, and numerical methods were developed to deal with the special mathematical 
features that would be common to all standing-wave, helically symmetric computations. These features 
include: (i) a mixed boundary value problem (regions of the domain in which the equations are hyperbolic 
and other regions in which they are elliptic) ; (ii) an iterative construction of nonlinear standing wave solutions; 
(iii) the extraction from the standing wave solution of an approximate outgoing wave solution. Reference (23j 
introduced a new technique that promised to reduce the computational burden: "adapted coordinates" that 
were well suited to the geometry of the problem both near the sources and far from them. In Ref. [2^ it 
was shown that with these coordinates good results could be computed by keeping only a very small number 
of multipoles, typically just the monopole and quadrupole moments. That method was applied to linearized 
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general relativity in Ref. (2j| and to the post-Minkowski approximation to general relativity in Ref. [25[. In 
the current paper we apply this method, PSW in adapted coordinates, to full general relativity. 

The mathematical infrastructure for this is rather involved. It includes the definition of standing waves 
for nonlinear fields; extraction of the outgoing solution; the adapted spatial coordinates; the multipole de- 
composition in adapted coordinates; the reduction to helical symmetry for tensorial fields; the full vacuum 
field equations; the inner boundary conditions on surfaces outside black holes. This infrastructure has been 
thoroughly documented in earlier papers, and the complete infrastructure is presented in Ref. [25j . Rather 
than give yet another lengthy presentation of this infrastructure, here we will give only enough of the math- 
ematical background necessary in order to explain the numerical results we are presenting, and to compare 
them with other numerical results. 

The rest of the paper is organized as follows. Sec. |TT] develops the mathematical description of helically 
symmetric tensor fields in terms of "helical scalars," quantities that are functions only of corotating coor- 
dinates. The field equations of general relativity, and of its weak field approximations are given for these 
helical scalars, along with boundary conditions on these fields. Section ITTTl discusses two important aspects of 
the numerical approach to our computations: (i) "discrete spherical harmonics" on our computational grid, 
and (ii) the use of the Maple^^ symbolic manipulation language to handle the complexity of the equations to 
be solved numerically. The numerical results produced by these methods are presented in Sec. IIVI and are 
discussed in Sec. El 



II. THE MATHEMATICAL PROBLEM 
A. The field equations and coordinates 

We start with the concept of a strong field source region and a weak field source region. The strong field 
source region is close to the inner computational boundary, surfaces with spherical topology around the black 
holes, and close to them but far enough from the horizon for computations to be feasible. We invoke a set 
of coordinates labelled t, x, y, z for which the spacetime metric 17^,^ deviates only slightly from the flat metric 
rjfj^i, at large coordinate distances from the binary sources. 

We impose helical symmetry by requiring that in these coordinates there be a Killing vector given by 

^^dt + n{xdy-yd^) (1) 

where 17 is a constant. It is convenient to define a rotating coordinate system by 

t = t z = z X = X cos fit + y sin fit y = —x sin fit + y cos fit . (2) 

and to introduce two cylindrical coordinates sytems: r, z, 4> in terms of x, y, z by the usual flat space formulas, 
and r, z, in terms of x, y, z also by the usual flat space formulas. We note that Eq. ([T]) is equivalent to 

^ = dt + fl {xdy - ydx) = dt+ fld^ = c^, (3) 

We follow the formulation of Landau and Lifschitz [2^ [S^l for the Einstein equations, which encodes the 
geometric information in the densitized metric 



g"'' ^ v^d^g"''. (4) 
To simplify the fleld equations, we choose to have our coordinates obey the harmonic condition 

5"^/3 = 0, (5) 
we define the inverse g^^^ of our basic field by 

S"'^g^^ = 5%, (6) 

and we define /i"^ by 



g"''^ y^drt^ -/i"'^). (7) 
The vacuum field equations then take the simple form 
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where 
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(25. 



(9) 



Equation ([H]), along with the definition ([7]), arc the equations to be solved for the unknown fields h""^ . 

Note that so far Eq. ([5]) is exact; there is no assumption that the Jf^ fields are perturbative. If we do treat 
these fields as perturbative, the post-Minkowski approximation follows by replacing g on the right by 77 to get 
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For linearized general relativity, the terms on the right can be ignord and the field equations are simply 
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These can be solved in closed form as expansions in special functions 24l|. 

In the case of a scalar field model, helical symmetry means that the scalar field $ is a function only of 
corotating coordinates x, y, z, or r, 0, z, and not of t. The computations for the field then involve only three 
coordinates. Complications arise when dealing with the components of tensor fields, since computational 
fields must be "helical scalars," i.e. , functions only of the rotating coordinates. The approach we employ is 
to define 
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in which [/("i)^i/("i), $(00)^ $(20)^ [7(21)^^(21)^ [^(22)^^(22)^ ^^.^ ^^^^^ helical scalars that 

carry all the information about the metric. 

We denote any of 4 real and 3 complex helical scalars as with A taking the value nn,nO,nl, .... In 
terms of these helical scalars the field equations ([8]) take the form 



- 2i^(A)f7^9y*^ + fi{Ayn^^'' = Q'^(*^) 



(19) 



where fi{A) has the value of for A = (nn), (nO), (00), (20), has the value ±1 for A= {n± 1), (2± 1) and has 
value ±2 for A = (2±2). 

In order to give the details of the source term we must introduce a set of objects n, e^., e^^, e^. These objects 
are described in detail in Ref. 25], and can roughly be considered to be the basis vectors associated with 
coordinates t,x,y, z. In terms of these, we define the objects 
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which are shown to be heHcal scalars in Ref. ^2^. With this notation the exphcit form of the source is [SC 



B ''C 



t ,f,C 



(27) 



Here the summation is not only over the tensorial indices, but also over the indices, B, C which range over 
the 10 values (nn), (nO), {n ± 1), (00), (20), {2 ± 1) J2 ± 2). 

Rather than work with the complex fields ^I'f"^), ij;(22)^ practice we work with the real and imaginary 

parts and V^, and for A = (nl), (21), (22) Eq. ^ is replaced by 



+ 2^{A)n'd^V + fiiAyn-'W = Real part of (Q^) 
nV^ - 2^i{A)Vi^d^U'^ + ii{Af^^V'^^ = Imaginary part of [Q^] . 



(28) 
(29) 



The □ operator here, as in Eq. (|19p. is ri°'^dadfj^ and the helically symmetric time derivatives are implemented 
through the replacement dt —il{xdy — ydx) = —ild^p, so that 



(30) 



The mathematical problem of finding the fields in full general relativity then consists of solving the partial 
differential equations (fT9| , (j28p , ((29|) . The post-Minkowskian and linear approximations to general relativity 
follow from making the appropriate simplifications of Q-^. 




FIG. 1: Two-center bipolar adapted coordinates. On the left is shown curves of coordinates x and O in the $ = 
orbital plane. On the right are surfaces of constant x, ©, and "!>. At large x/'^ the adapted coordinates 0,<I> are 
spherical polar angles with respect to the corotating coordinates X, Y , Z, a permutation of the coordinates x, y, z, 
discussed in the text. 

Computational solutions of the field equations are carried out in an "adapted" corotating coordinate system 
X,9,$ illustrated in Fig. [TJ These coordinates are two-center bipolar coordinates, defined in terms of the 
corotating coordinates x, y, z. The two "centers" are the coordinate points x = ±a, y = z = 0, points that 
roughly represent the locations of the black hole sources. At large values of x/oj the adapted coordinates 
approach ordinary spherical coordinates with x the radial coordinate. At small values of x/a, the adapted 
coordinates become a reparamctrization of spherical coordinates centered on the bipolar points, with the 
X coordinate approaching y2aJZ where TZ is the (corotating) coordinate distance from one of the bipolar 
centers, and with 20 approaching the angle with respect to the x axis. These coordinates are discussed in 
detail in Ref. [2^. 

Each of the helical scalar fields is expanded in the form 

* = E '^'i^ix)YUe,'^)- (31) 

em 

Due to the nature of the coordinates only a few multipoles need to be kept. The argument for this "multipole 
filtering" is given in Ref. [11] , along with numerical results demonstrating its validity. 
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B. The boundary conditions 



The standing wave condition is a procedure rather than a constraint, hke Dirichlet, Neuman or Sommerfeld 
conditions. This procedure is apphed in the weak wave zone at some value Xmax of x that is many times 
larger than a, and much larger than a wavelength (1/17). In this weak wave zone the fields can be treated 
as perturbations of flat spacetime, and outgoing/ingoing conditions on the perturbations can be imposed. 
The standing wave procedure for full general relativity starts with the system of equations (|1]) ~ ([5]). An 
approximate solution for h"^ is put into the right hand side of Eq. ([8]), both for the explicit appearance of 
the /i"^ terms and through the dependence of the g"^ terms on h"^ . The right hand side is then treated as 
known, and Eq. ([8]) is considered to be a linear equation for /i"^. In Ref. [23] it was shown that in adapted 
coordinates outgoing (+) and ingoing (-) conditions on any of the helical scalars are 

1 9 / ~\ ^ ( cose . 

-7^ X* = cos*- s n$— . 32) 

X9x V-^ y \^ 99 sine j ^ ' 

These conditions are used to compute both outgoing and ingoing solutions to the linear equations for each 
of the helical scalars Vl/. The ingoing and outgoing solutions are then averaged, the result is taken as an 
improved approximation for the field, and the process is iterated. 

The convergent result, which we call the standing wave solution, is an exact (numerical) solution of the 
nonlinear field equations with no net flow of energy inward or outward. The field equations could also be 
solved for a nonlinear outgoing solution, or a nonlinear ingoing solution. In principle, for nonlinear equations, 
the standing wave solution is not the average of the ingoing and outgoing solutions. In practice, however, it 
is an excellent approximation for the following reason, which we call "effective linearity." The ingoing and 
outgoing solution in the strong field region is negligibly affected by the boundary conditions in the strong 
field zone, so that the average of the ingoing and outgoing solution there must be the same as the standing 
wave solution. In the weak wave zone, the solution is completely different for different radiation conditions, 
but in the weak wave zone the problem is effectively linear, and again the average and the standing wave 
solutions are very nearly the same. 

Effective linearity means that we can treat the standing wave solution as if it were the average, and we 
can decompose it into approximate ingoing and outgoing solutions. We take the outgoing solution to be our 
approximation to the solution to the physical problem. 

The inner boundary is imposed on quasispherical surfaces, around the bipolar centers, at some small value 
Xmin of X- The choice of the numerical value of Xmin is a compromise. The choice of a small value means 
that the conditions at Xmin will be dominated by the presence of the black hole it surrounds, and the other 
hole in the binary pair can be considered to have a weak infiuence. But a small value of Xmin means that the 
computational regime Xmin < X ^ Xmax , includes fields strong enough to cause problems in the convergence 
of the computational method. 

An important point about the inner boundary conditions is made in all previous papers on this method: 
the details of the inner boundary conditions are equivalent to the details of the sources they represent inside 
the X — Xmin surface. But the details have almost no effect on the fields in the wave zone, and even in the 
intermediate x ~ ^ zone. We can, therefore, impose conditions at x = Xmin that represent some reasonable 
source, and defer the task of adjusting the details. The wave fields, the computed outgoing energy, etc. will 
be negligibly affected. 

Our choice is taken to be the conditions for a single Schwarzschild hole comoving with harmonic coordi- 
nates. This reasonable choice was used in our post-Minkowksi computations, and hence that choice facilitates 
comparisons with previous work. The actual form of the boundary conditions have been previously reported 
[25| and are included here for completeness: 

/4M 7M2\ 2 _ M^v^2<^ ^^^2 _2 . 



^^t.-^(2e)cos^'J> (33) 



_ /l,/-2 2 4 

vI/"° = ± — ^|^sin2(2e)sin<i>cos$ (34) 



U^"'^ = |J^^«in(2e)cos2ecos$ (35) 
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4M TAP 



sin2(2e)cos^$ 



sgn[cos 9] 



(36) 



(37) 



^20 ^ 



(cos^ 29 + 7^ sin^ 29 cos^ $ - 2 sin^ 29 sin^ $) (38) 



U 



7^4 4a2 



sin 29 cos 29 sin $ sgn[cos 9] 



V 



7^4 4a2 



sin 29 sin $ cos $ 



[7(22) = _ rii!i 



4A/ 7M2\ t;2^2 ^2 ^4 



7^ 7^2 / 2 27^4 4a2 



(cos^ 29 - 7^ sin^ 29 cos^ $) 



(39) 
(40) 
(41) 



y (22) = _ COS 29 sin 29 cos $ sgn[cos 9] 



Here, TZ, in terms of adapted coordinates, is given by 



(42) 



4a2 



[1 +72i;2sin2 29cos2$] , 



(43) 



where v = afl is the coordinate speed of the sources and 7 is the corresponding Lorentz factor — f2 . 

The symbol TZ is the coordinate distance from a source "point" in a frame comoving with the source; its form 
in adapted coordinates is derived in Ref . [2^ . The ± in Eq. ([M]) distinguishes betweeen the condition (+) for 
the source point at 9 = and the condition (-) for the source point at 9 = tt. 



C. Numerical implementation 

In practice one of the issues that requires considerable care in numerically implementing the PSW method 
is treatment of the boundary conditions. Consider first the outer boundary conditions. For the most part, 
numerical implementation of boundary conditions p2p consistent with expansion (j3ip is straightforward. 
However, special care must be taken with the non-r adiative parts of *(20)^ ^^^^ ^(oo) ^ ^^^^i the 

analytical and numerical solutions of these fields are dominated by non-radiative modes. In Ref. [l^l the 
general form of the analytical solution for thes fields in linearized theory is reported to be 



2i+l' 
e=o,2,... 



iXn Y mje{mnr<)Y;^{Tr/2,0)Ye„^{e,0)lm\h[^\mnry)e'"''^] . (44) 



=2,4,... m=2,4,. 



Here ii' is a constant that depends on whether A is (nn), (20), or (00); ji{x) is the spherical Bessel function 
of order £; ft,^^-* is the spherical Hankel function of order £; r< = min(r, a), and r> = max(r, a). 
Expression shows that series solutions for ^("^"), ^(2'0)^ and ^^o^") 

are dominated by the £ = 0, m = 

mode, which decays as 1/r for large r values. Since the radiative part also decays as 1/r it is necessary to 
separate numerically the radiative and non-radiative parts of the solution at large radius. This is accomplished 



7 



by setting the following boundary condition for the monopole apo (i-e., the £ = 0, m = mode amplitude) of 
these fields 

— aoo = ■ (45) 

dx X 

This is just a mathematical statement that at large distances the dominant part of aoo falls off as 1/r. 

The numerical implementation of the inner boundary conditions must also be treated with care. It turns 
out that some of the conditions in Eqs. ([55]) - are not to be considered as explicit values to be set at the 
inner boundary, but rather as regularity conditions on certain of the fields. (It should be understood that 
not all the fields need to have data imposed to give them a scale. Due to the coupling of fields in the field 
equations some of the fields are given scale by the coupling to other fields, not by their own boundary data. 
In practice we do this by replacing Eqs. (^0]) . and (H^ with 





= 


(46) 


dx 


= 


(47) 


dx 


= 


(48) 


Ay(22) 
dx 


= 0, 


(49) 



since this choice has already been used in Ref. [2J| and [2£ 



III. NUMERICAL METHOD 
A. Discrete spherical harmonics 

Discrete spherical harmonics were introduced first in Ref. [24] as an integral part of the eigenspectral 
method. In a straightforward approach to their use, spherical harmonics of analytic theory would be evaluated 
at discrete grid points. We have found that this leads to unacceptable failures of orthogonality, so instead we 
have used "discrete spherical harmonics" which are orthogonal to the level of machine precision. 

To understand these, consider a two dimensional computational grid with N = hq x n$ points. The discrete 

(k) 

spherical harmonics Y^j are A^-dimensional vectors which are conveniently represented by the two indices 
1 < j < riQ and 1 < j < n,^. Let the matrix Lab,ij be the operator equivalent to the angular Laplace operator 
evaluated at 8a and $h . This operator is therefore defined by 

[sin evLg/(e, $)] « ^ L„fc,.,/,, (50) 

where /(0,$) is an arbitrary function and where fij — f{Qi,^j). The approximation symbol in the rela- 
tionship is due to the truncation error induced by the finite difference representation of the derivatives on the 
left. 

(k) 

The discrete spherical harmonics Y-j ^ are defined to be the solutions of the generalized eigenvector problem 

E Lab,^JY}'^ - -A Sin BaV!^^^ , (51) 

together with the normalization condition 



EE^f^f ^0^*-^^^'- (52) 

i 3 

Details of the construction of the matrix operator Lab,ij and of the computation and properties of the eigen- 
vectors can be found in Refs/ [13] and [soj . 
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B. The use of symbolic manipulation 

It is clear that most of the aspects involved in the computation of fully relativistic fields with with the 
eigenspectral method (i.e., with adapted coordinates, multipole filtering, and discrete spherical harmonics) 
are not conceptually difficult, but are technically difficult to implement. The advantages of the eigenspectral 
method, explained in detail in references [H, [IBi offset by an increased complexity in the form of the 
equations. The several stages required by a typical algorithmic solution to a problem within the eigenspectral 



method are illustrated in Fig. 2(a) During the early development of the PSW program most of the work at 
these stages was carried out by a human being, with the computer used only for the final numerical solution 
of the equations. 



The fact that most of the stages depicted in Fig. 2(a) are algorithmic in nature suggested the idea of 



developing a computational framework to deal with the different aspects of the computation. The framework 
was designed to ease the numerical implementation of any level of approximation within the PSW model 



(i.e., linear, PM or fully relativistic) by reflecting the workflow of Fig. 2(a) This framework included the 
implementation of several sets of tools for the PSW project within Maple ^ , a general-purpose computer 
algebra system (CAS) software package. 

Three different set of tools were developed for the PSW project. The first set corresponds to the conversion 
of an equation in Minkowski coordinates t, x, y, z into adapted coordinates. The role of a second set of tools 
was to perform the conversion to finite difference form of any given differential equation. Finally a third 
set corresponds to tools implemented to convert Maple expressions into C functions. This final set includes 
tools that are able to create all the different C functions related to the projections over spherical harmonics 
and the construction of the matrix system that has to be solved at each iteration either in a perturbative 
or exact scheme as described in Ref. [1^. The code generated by these tools was later embedded into a 
bigger infrastructure developed in C. The C infrastructure essentially takes care of runtime issues, such as 
the allocation of memory and the interface of advanced numerical routines in LAPACK for matrix inversion 
and eigenvector computation for the routines generated by Maple. 

The implementation process for a given model within the PSW-eigenspectral method (namely, linear gravity, 
post-Minkowski, or full general relativity) is streamlined through these tools. First the differential equations 
have to be provided in completely explicit form. The first set of tools is then applied to these equations, 
rendering the equations in adapted coordinates. The equations that are output are next passed through the 
second set of tools, to put them into finite difference form. Finally, the output of this process is passed 
through the third set of tools in order to obtain the C code. 

The net effect of the use of these tools "raises the bar" for the implementation of numerical models with the 



eigenspectral method, in that most of the process is performed by the computer, as Fig. 2(b) illustrates. We 
emphasize that the use of symbolic manipulation was indispensable to the last phase of the PSW-eigenspectral 
program. This leads us to suggest that a similar eclectic use of algebraic manipulation software and numerical 
programming is a powerful approach not only for our problem but for other problems in which code must 
be generated for complex and lengthy mathematical expressions for which human coding is prone to error. 
However it is worth noting that neither Maple nor any other computer algebra system is able by itself 
can directly generate the C expressions required for this kind of project. Although most algebra systems 
features parsers to generate C code most of the time the native parsing algorithm can not deal with some 
subtleties. Some of these subtleties include type specification for intermediate variables, data type signature 
of functions, the appropriate translation of some mathematical functions and the creation of header files 
among others. For the present work it became necessary to modify the native parser of Maple in order to deal 
with these subtleties. In particular the syntax definition of the C parser was enlarged in order to include some 
mathematical functions while the parsing algorithm was modified in order to ensure the correct signature for 
the functions defined during the automated process. In addition some fine tuning was performed in the parser 
algorithm in order to avoid incorrect or unnecessary operations in the translated code, such as some type 
casts or incorrect array indexing. These kind of subtleties must be addressed and tailored to the occasion by 
any other project that wants to use a similar approach of the mixing Computer Algebra Systems and analytic 
tools. 



IV. NUMERICAL RESULTS 

The ultimate unknowns for which we actually solve are the coefficients a^j(x) of expansion of the 
fields in Eq. (PT|) . The equations for these unknown a^(x) are the field equations (fTO|) . along with 
the appropriate inner and outer boundary conditions. In the source term for the field equations, the 
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FIG. 2: The process of implementing numerical modeling has been modified by the use of Maple. The aid of a computer 
algebra system has allowed the computer to handle a great deal of the "analytic" work, significantly simplifying the 
role of the programmer in the development of numerical code 



dependence on af^{x) occurs both in the explicit quadratic appearance of ^I''^ on the right in Eq. (|27p . and 
through the very compHcated dependence of S"^^^ on encoded in the equations ([7]) and p^ - P^ . 
In our work on the post-Minkowski approximation [25j only the first dependence was kept, since inclusion of 
the dependence through 5*"^^^ would have gone beyond the order of the approximation. But in the present 
work, for the full equations of general relativity, we cannot ignore this deeper nonlincarity. In fact, it is 
precisely this dependence that, in our formalism, distinguishes full general relativity from the second-order 
post-Minkowski approximation. The complexity of this nonlinearity, however, leads to a practical problem. 

In our work, here and previously, with nonlinear equations, the mathematical problem to be solved was 
cast into the form L{y) = F{y) in which L is a linear operator on the unknown function, or set of functions y. 
A solution could be sought in two ways. A direct iteration could be used in which an approximate solution 
yn is substituted into F(y) and the linear equation L{yn+i) ~ F{yn) for a new approximation yn+i is solved, 
with appropriate boundary conditions. Alternatively, Newton-Raphson iteration can be used. In this method 
a linear equation is found by expanding the y dependence of F{y) about _F(y„), with a Jacobian playing the 
role of the derivative of F with respect to the fields y. Newton-Raphson iteration is expected to have better 
convergence properties than direct iteration, and this was indeed found to be the case in previous work. 

In principle, then, it is a Newton-Raphson scheme that is to be sought. But due to the complicated depen- 
dence of S"^^^ on a^(x) this does not turn out to be feasible. With Maple generated code the evaluation 
of the Jacobian needed in the Newton-Raphson method becomes unreasonable long. The main reason for 
this is that the code generated by Maple is optimized for execution time by unrolling lengthy evaluations 
explicitly in several steps using several temporary variables. As a consequence this highly optimized code is 
extremely lengthy in terms of the number of code lines. As the number of terms given to Maple increases so 
does the code that must be compiled. If the number of terms is high enough the raw code generated with the 
Maple tools becomes so lengthy that it is impossible to compile. This has severely limited the applicability of 
this automated process for Newton-Raphson methods for a fully relativistic model. For this reason the fully 
relativistic model was ultimately implemented using a direct iteration scheme. 

We have found that the use of direct iteration does not seem to be a major problem. In particular, 
models with parameters that allowed convergence for post-Minkowski models with Newton-Raphson iteration 
have been found to converge in full general relativity with direct iteration. One reason for this, is that 
we have limited our application of the full general relativity models to those for which the post-Minkowski 
approximation is valid. As explained in Ref. (25j . this means that Xmin must not be chosen very small, to 
avoid very strong nonlinearities. This in turn leads to the constraint on model parameters [2^ 

2V2an < Xmi„/a < V2 . (53) 

With this restriction on field strengths convergence has been found to be very good. 

This successful convergence is illustrated in Fig.[3]for ^'(2:0)_ These results show that convergence is attained 
after 3 to 4 iterations for angular velocity Q = 0.15; similar or faster convergence is found for all smaller values 
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FIG. 3: Typical convergence of fully relativistic using direct iteration. The field is shown as a function of 

coordinate radius along a line in the orbital plane approximately through the source points (i.e., at O = AO). The 
simulation parameters here are = 0.15, x ne x n<j = 1750 x 16 x 32, with Xmin/o = 0.50. The computation 
includes all multipoles up to ^ = 3 (monopole and quadrupole terms kept) although the monopole term aoo is not 
included in the results displayed. 

of Vl. We usually consider = 0.15 to be the limit of the angular velocity for which the PSW method gives 
a good approximation to the outgoing solution. 

A very important consideration for the success of the computations, and of the PSW approach, is the 
value of Xminj since this parameter determines the maximum field strength encountered in the computation. 
Moreover, results that are sensitive to the choice of Xmin raise questions about the underlying assumption in 
the PSW approach that the details of the inner boundary condition are not important to the fields except 
very near the inner boundary. Figures H] show the results for all fields, for a model with — 0.15. 

These results, show somewhat mixed results for the insensitivity to the value of Xmin- For the "coulombic" 
field ^f^"'") the computed result is particularly insensitive to Xmin- Several of the fields show significant 
diffrences in the results for the three values of Xmin- One would expect that the results for the largest value 
of Xmin should havc errors in representing the boundary conditions appropriate for a "point" particle, since 
Xmin/a = 0.6 is not really justified as a "near source" surface. One might expect better agreement between 
the two smaller values of Xmin, since both inner boundaries are reasonably close to the source "points." The 
results tend to confirm this, though with exceptions. Of particular interest are the results for [/(2,2)^ y(2,2)^ 
since these fields carry the information about gravitational waves in the weak wave zone. For these fields the 
results for Xmin/a = 0.45 and Xmin/o = 0.3 agree moderately well in phase, but have amplitude differences 
as large as a factor of 2. A study of the origin of this difference will help clarify the nature of the correct 
boundary conditions. It should be noted that the computations of the fields within full general relativity have 
turned out to be less sensitive to the choice of Xmin/a than are the post-Minkowski computations. In some 
way, the inclusion of the dependence of iS"^!^^ on the fields moderates the sensitivity of the results to Xmin/o- 

Another type of comparison is presented in Figs. [5l the comparison of computations based on linearized, 
post-Minkowski and fully relativistic models for O = 0.075 and Xmin/a = 0.4. The monopole term aoo is not 
included in the results shown for ^fCO'"), \|f(2,o)^ g^j^^ The results show that and are very 

well approximated within the post-Minkowski model. All other fields differ in more or less a significant way 
with respect to the post-Minkowski results, showing that the fully relativistic contributions are of significance, 
even for the moderate field strengths of these computations. In this connection it should be noted that the 
field is very well approximated by lower order approximations. It is only this field that was used in 

the source term for the iterations in Ref. [l^. Thus, it is not the imperfections in that are the basis 

of the inaccuracies, but rather the suppressed dependences on additional fields. 
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FIG. 4: Comparison of results for direct iteration of relativistic computations as function of Xmin for SI — 0.15. Fields 
are shown in the orbital plane along a line approximately through the source points. These results reveal that the 
location of the inner boundary is less important in fully relativistic models than in post-Minkowski models described 
in [H, specially for 
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FIG. 5: Comparison between linearized, post-Minkowski and fully relativistic computations with Q — 0.075. The 
computations used a numerical grid of x ne x n<s> = 1750 x 16 x 32, with Xmin = 0.40a and Xmax ~ 100a. 
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V. DISCUSSION AND CONCLUSIONS 

The results presented in the previous section for the fuUy relativistic case show the same convergence rate 
as was found for the directly iterated post-Minkowski approximation and show similar limitations in field 
strength as the Newton-Raphson method for the post-Minkowski approximation. The results also reveal 
that the nonlinearities present in the fully relativistic model produce important changes in the numerical 
solutions for most of the fields. For the most important field of the post-Minkowski approximation ^("'"), 
it is comforting^ to note that there is excellent agreement of the post-Minkowski and fully relativistic results. 
The fact that is also well represented by the post-Minkowski computation is also a good sign. These 

good agreements are expected qualitatively since the post-Minkowski approximation is correct up to second 
order in the small parameter M/TZ while v[((".") and vP^"'^) are of order 1 and 1.5, respectively, in the same 
parameter. In this sense the relativistic results confirms the validity of the post-Minkowski approximation. 

It is important to mention not only the strengths of the method but also its shortcomings. First, the use 
of adapted coordinates gives a solution on an irregular computational grid on the Minkowski background of 
the problem. As a consequence the results obtained from the computations require a difficult interpolation 
if they are to be compared with other results or used as initial data for evolution codes. Second, the method 
is inherently limited in its accuracy. Higher accuracy, would ultimately require more multipoles, and an 
increase in the number of multipoles adds to the size of the equation set to be generated by Maple, and to 
the compilation problem. 

The current accuracy, however, is sufficient for us to see in the results presented in the previous section that 
there is sensitivity of the computed fields to the location of the inner boundary Xmax- The inner boundary 
conditions used for the computations are those for a a spherically symmetric mass point, or Schwarzschild 
hole, boosted so that it is instantaneously comoving with the rotating coordinates. (For a detailed discussion 
see Ref. [HI). This is, of course, an approximation that ignores the tidal effects of the two source points on 
each other and the fact that the sources are not moving in straight lines of a Minkowksi background. That 
approximation would be adequate at sufficiently small il, but - the results suggest - not at fl = 0.15/a. 

Improved inner boundary conditions are an important application of the PSW method, since this is part 
of the use of the method to study radiation reaction. But the issue of the inner boundary condition may not 
be separable from the issue of the gauge condition, which might explain the sensitivity of the results to Xmax 
even for the relatively small angular rotation rate fl = 0.15/a. The "gauge issue" is the fact that we have 
used a specific gauge, that of Eq. ([S]), but we have not enforced or checked that gauge in the solution. If the 
problem is well posed then there should be a solution which can be transformed to this gauge. The question 
of well-posedness directs attention to the inner boundary condition, and the multipoles at the inner boundary. 
The coarsest parameter in those boundary conditions is simply the mass that we ascribe to the source, a mass 
that we have linked to the angular velocity through the nonrelativistic Kepler law. One would certainly 
expect relativistic corrections to that relationship, so in a sense our solution of the fully general relativistic 
problem has omitted an important source of nonlinearity. 

With the adapted coordinate method we have demonstrated that the PSW computations can be brought 
to a level of accuracy adequate for studying features of the strong fields and inner boundary conditions. It 
is therefore capable, in principle, of giving useful insights about radiation reaction contained in the details 
of the near fields. The method can also be used in the future to provide a quasistationary sequence of 
simulation parameters (i.e., mass, orbital radius and angular velocity). Such sequence will represent more 
than adequately the slow inspiral process of the binary. The main challenge for this approach is that it 
is important to show that we are comparing the "same" system at different radii (i.e., that some physical 
quantities are conserved although energy is lost during the process). For binary neutron stars the baryon 
number must be unchanged. In the case of black holes the issue is more difficult, since the total mass of the 
system decreases as energy is emitted. However the local mass of the black holes should not change; this local 
mass could be calculated at each stage of the sequence using numerical codes designed to detect horizons 
locally, as described by Ashtekar et al [1^. 

The current work has solved the problem with the accuracy that was sought, and a broader lesson might 
lie in that success. The success was achieved through the combination of a Computer Algebra System such 
as Maple and numerical programming. This combination can be a powerful approach to some problems that 
requires not only fast and sophisticated numerical algorithms but also a fast and sophisticated implementation 
process. The lesson in the PSW program is that the use of a highly analytical approach to the problem 
(adapted coordinates and multipole filtering) translated into an increased complexity in the algebraic form of 
the equation. While humans are prone to errors when manipulating lengthy expressions, computers are not. 
Leaving the details of the manipulation and coding to an algebraic manipulation language is a natural choice 
for this sort of problems. Of course such an approach involves the development of tools adapted to the needs 
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of a particular project within Maple or any similar development platform, which can be quite challenging. 
However, once the tools are developed and benchmarked the correctness of the numerical implementation is 
guaranteed, which is a highly desirable feature for projects involving complex numerical implementations. 
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